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I.  IKTROUCTION; 


Attempting  to  determine  fluid  flow  around  an  obstacle  requires  the 
solution  of  the  Naviex-Stokes  equations; 


where 
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Here,  the  four  dependent  variables  p,  u,  v,  and  e  represent  the 
physical  quantitive  of  density,  x-  and  y-  components  of  velocity,  and 
Internal  energy. 

Even  a  casual  glance  at  this  system  of  equations  shows  that  It  has 
many  extremely  difficult  aspects.  It  Is  a  system  of  four  equations  in 
the  four  basic  unknowns,  which  Is  of  mixed  parabolic  hyperbolic  type. 

It  Is  In  two  space  dimensions  and  In  addition  Is  nonlinear.  Any  one  of 
these  problems  would  be  sufficient  to  make  obtaining  an  analytic  solution 
difficult  so  It  Is  fair  to  say  that  an  analytic  solution  Is  out  of  the 

t 

question.  Thus  these  problems  are  "solved"  by  numerical  methods. 


What  this  means  of  course  is  subject  to  some  interpretation.  What 
usually  happens  Is  that  a  nuoiber  of  approximations  are  made  to  the  original 
problem,  in  order  to  arrive  at  a  discrete  approximation  which  may  then  be 
solved  by  computer.  Most  obviously,  'space''  Is  Interpreted,  not  as  a 
continuum,  but  as  a  finite  number  of  grid  points.  The  value  of  pressure 
for  example  at  a  grid  point  Is  then  used  to  calculate  the  value  at  nearby 
points,  by  means  of  a  few  terms  In  a  Taylor  series.  Derivatives  are 
approximated  by  finite  differences.  If  the  flow  Is  around  a  body,  then 
one  ass\imes  that  velocity  is  zero  on  the  surface  of  the  cylinder,  and 
the  surface  temperature  Is  prescribed.  The  assumption  Is  that  If 
the  grid  Is  fine  enough,  and  If  time  steps  are  sufficiently  small,  then 
this  discrete  model  Is  a  reasonable  approximation  of  the  original  eq¬ 
uations  (1) . 

A  method  for  solving  supersonic  flows  that  was  found  by  experience 
to  work  well  was  MacCormack's  alternating  direction  explicit  scheme. 

Thus,  this  was  a  logical  method  to  be  experimented  with  for  subsonic  flow. 
However,  the  transition  was  not  easy.  The  MacCormack  method  demands 
certain  fictitious  (or  numerical)  boundary  conditions  due  to  the  dif¬ 
ference  algorithm  which  are  not  physically  present,  and  which  were 
arrived  at  by  a  process  of  computational  experimentation.  These 
methods  did  not  work  In  the  case  of  subsonic  flow,  although  they  did  In 
the  case  of  supersonic  flow. 

II.  OBJECTIVES; 

The  purpose  of  this  project  was  to  explain  thus  anomaly,  and  to 
arrive  at  Improved  far-fleld  boundary  conditions  for  the  subsonic  case. 

To  do  this,  we  will  first  review  the  theory  for  a  single  wave  equation. 
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then  for  systems  of  linear  hyperbolic  partial  differential  equations,,  and 
finally,  we  shall  discuss  how  numerical  experiments  with  the  Navler-Stokes 
equation  confirm  the  predictions  based  on  the  elementary  theory. 


Ill  REVIEW  OF  NUMERICAL  THEORY  FOR  ONE  EQUATION; 

Obviously  the  full  system  of  equations  (1)  Is  too  complicated  to 
achieve  a  great  deal  with  sophisticated  mathematical  analysis.  Typically, 
mathematicians  will  work  with  simpler  equations  which  In  one  way  or 
another  resemble  the  original  system.  One  tries  to  use  the  mathematical 
Insights  of  the  simpler  situation  In  the  context  of  the  more  complicated 
one.  This  is  obviously  not  one  unbroken  chain  of  reasoning,  but  more  a 
process  of  educated  guesswork. 

Surprisingly,  one  equation  which  provides  a  great  deal  of  Insight 
Is  the  simple  wave  equation 


u  +  au  -  0 


(2) 


We  consider  this  equation  on  the  Interval  [0,ll.  On  the  real  line,  the 
solutions  to  this  equation  would  be  waves  running  from  left  to  right  with 
velocity  a,  which  are  constant  along  the  lines  x  -  at  *  d,  deR 

Thus,  If  we  consider  this  equation  on  the  region  {(x,  t)  0^  x  ^1,  t^O} 
then  It  will  be  analytically  determined  by  the  Initial  value  U  (x,  0)  >  f(x 
O^x  <^1  and  the  boundary  value  u  (o,  t)  •  g(t).  With  these  conditions, 
the  Initial  boundary  value  problem  Is  well  posed.  It  *ls  Impossible  to 
prescribe  boundary  conditions  at  x  *  1  Instead  of  x  >  0  without  either  (a) 
limiting  the  Initial  values  or  (b)  causing  discontinuities.  All  of  this 
makes  perfect  physical  sense.  In  order  to  know  what  happens  In  a  wave 
sltuatlCT),  we  must  give  Information  on  the  Initial  conditions  and  also  on 
the  waves  entering  at  the  Inflow  point. 


However,  the  MacCormack  scheme  which  Is  equivalent  to  the  Lax-Wendro£ 
scheme  requires  some  knowledge  of  u  (1,  t)  Indeed  the  value  for  Uj^^  is 
given  in  terms  of  Uj_2,  Uj*  Therefore,  the  boundary  point  x  -  1, 

(J  -  J)  requires  special  treatment  of  some  sort.  This  special  treatment  is 
called  a  "numerical  boundary  condition*'  or  compatibility  condition.  The 
imposed  can  have  an  enormous  Impact  on  the  successful  numerical  solutions 
of  the  problem. 

Pherhaps  the  best  two  recent  summaries  on  this  simple  equation  were 
accomplished  by  Kreiss  [s]  and  Gottlieb  and  Turkel  [2]. 

Krelss  is  basically  concerned  with  what  can  go  wrong.  He  observes 
for  example  that  if  you  overspecify,  l.e.  just  make  a  guess  at  what  u(l,  t) 
is  going  to  be  and  then  simply  prescribe  it,  convergence  to  a  steady  state 
may  or  may  not  take  place,  depending  on  whether  the  number  of  grid  points 
is  even  or  odd.  If  this  happens  with  this  simple  case  obviously  you  don't 
want  to  try  it  in  a  more  complicated  case.  (We  shall  have  more  to  say  on 
this  later). 

One  method  that  works  well  is  to  use  u^  (1,  t)  -  0,  or  in  finite  dif¬ 
ference  form  Uj  >*  Uj_^.  A  slight  error  in  u  at  the  outflow  point  is  made, 
but  since  the  flow  is  from  left  to  right,  this  error  does  not  propagate  back 
into  the  x-domaln.  This  is  proved  analytically  in  Parter  [S].  About  the 
worst  mistake  that  can  be  made  is  not  to  specify  u  at  the  incoming 
boundary.  For  example,  one  might  confuse  the  Inflow  and  outflow  boundary 
and  prescribe  u"  ■>  u^,  and  Uj  «  M  all  n.  In  this  case  as  the  space  step 
and  time  step  become  small  this  conveys  to  a  steady  state  on  an  Interval 
[0,T]  where  the  steady  state  is  determined  by  the  initial  conditions  at 


the  inflow  point. 


A  moment's  reflection  ought  to  convince  us  of  how  undesirable 
this  is,  since  we  do  not  know  the  correct  steady  state,  our  initial 
conditions  are  bound  to  be  different  from  the  correct  solution.  However, 
what  we  seem  to  converge  to  actually  depends  on  the  choice  of  the  initial 
conditions.  Again,  we  shall  have  more  to  say  about  this  when  we  discuss 
computational  solutions  of  the  Navler-Stokes  equations. 

Another  excellent  paper  on  the  same  subject  is  by  Gottlieb  and 
Turkel  C2].  In  this  paper,  a  comprehensive  review  of  many  different  num¬ 
erical  boundary  conditions  is  given.  For  example,  in  addition  to  the  ones 
already  mentioned  at  the  outflow,  one  might  consider  u^^  ■  0,  which  numerically 
is  u^  -  -  2u^  .  +  u^“  Oor  even  u  +  au  “0  where  »  u^  a  (At/Ax) 

(Uj  ^  -  Uj) .  This  corresponds  to  a  one-sided  ''upwind"  difference  approxi¬ 
mation  at  the  outflow  boundary. 

Their  conclusion  is  that  the  upwind  difference  appears  best,  although 
the  convergence  of  the  scheme  with  =  0  is  just  as  fast  (but  less  accurate). 

IV.  WELL-POSED  BOUNDARY  CONDITIONS  FOR  LINEAR  SYSTEMS; 

The  simplest  systems  of  two  linear  wave  equations  to  consider  is 
u^  +  au  »  0 

t  X 

V^  -  bv^  -  0  (3) 

where  a^^  >  b  >  0,  0  ^  x  ^  1,  and  t  ^  0.  In  this  case  the  waves  in  u 
travel  from  left  to  right  with  velocity  a,  and  the  waves  in  v  travel  from 
right  to  left  with  velocity  b.  Clearly  u  must  be  given  initially  and  at 
the  left  boundary;  whereas  the  values  of  v  must  be  given  at  the  right  bound¬ 
ary.  This  is  analytically  necessary  in  order  to  have  enough  information  to 
solve  the  problem.  However,  there  is  a  minor  complication.  The  incoming 


value  of  u  may  be  given  at  x  =  0  in  terms  of  v  (which,  after  all, 
is  determined  there) ,  and  the  incoming  value  of  v  at  1  may  be  given 
in  terms  of  u.  Thus,  the  initial  boundary  value  problem  (3)  is  well 
posed  if  the  initial  conditions  are 

u  (x,  0)  -  f(x)  V  (x,  0)  -  g(x) 
and  the  boundary  conditions  are 

u  (o,  t)  =  F^(t)  +  Cj^  v(0,  t) 

V  (1,  t)  *  G^(t)  +  u(l,  t)  (4) 

Notice  that  if  c^^  =  c^  ■  0,  then  the  problems  are  completely 

uncoupled,  and  each  is  a  copy  of  the  single  equation  in  section  IV. 

For  a  moment,  consider  the  case  where  Fj^(t)  *  Gj^(t)  “  0,  and  c^^  “  C2  *  1* 

In  this  case,  we  get  a  different  phenomenon.  Let  us  assume  that  u 

is  identically  zero  initially,  and  that  v  is  identically  1.  Then  the 

square  wave  in  v  travels  to  the  left,  causing  u  to  be  non-zero  at  the 

in-flow  boundary,  and  thereby  propagating  from  left  to  right  with  velocity 

1.  When  u  reached  x  «  1  it  would  in  turn  Influence  v  by  the  boundary 

condition  (4)  and  a  new  disturbance  would  propogate  in  v  from  right  to 

left.  This  is  the  reason  that  this  boundary  condition  is  called  a 

"reflective  boundary  condition” . 

Notice  that  if  both  u  and  v  travel  from  left  to  right,  i.e.  if 

u^  +  au  “0 
t  X 

v^  +  bv  »  0 

t  X 

where  a  >  0,  b  >  0,  0<^x_<l  and  t  ^  0,  then  this  problem  does  not 
arise.  The  initial  conditions  of  u  and  v  must  be  specified  as  must 
the  boundary  conditions  at  x  ■  0.  No  reflections  or  coupling  is  allowed 
to  take  place. 
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(5) 


For  a  general  hyperbolic  system  of  the  form 
+  Au  *0 

t  X 

one  must  diagonalize  the  matrix  A,  l.e.  one  must  find  a  matrix  T  so  that 
T  A  -  D 

where  D  Is  a  diagonal  matrix.  One  then  makes  the  substitution  W  ■  TU  and 
the  equation  (5)  transforms  to 
+  D  W  -  0 

t  X 

The  number  of  positive  eigenvalues  of  D  Identifies  the  right-running 

I 

variables  and  the  negative  ones  Identify  left  running  variables.  At  this 
.point,  we  emphasize  that  the  only  way  to  be  familiar  with  the  wave  nature 
of  (5)  Is  to  look,  not  at  the  physical  variables  u  In  which  the  problem 
was  originally  presented,  but  Instead  at  the  new  variables  W  which  are 
linear  combinations  of  the  old  ones.  Only  then  Is  the  wave  structure 
apparent . 

If  k  Is  the  number  of  positive  eigenvalues,  and  If  represents  the 
k-vector  of  these  k  coordinates  of  Wj  and  Is  the  n-k  other  coordinates 

then  unless  and  are  Independently  prescribed  at  the  left  and  right 
boundaries  respectively  then  we  may  get  reflective  boundary  conditions. 

For  well-posedness ,  It  Is  sufficient  that  at  x  >  0,  be  given  (possibly 
depending  on  and  that  at  x  >  1,  be  given  with  possible  dependence 

on  W^.  Only  then  can  acceptable  boundary  conditions  be  formulated  In  terms 
of  the  physical  variables  U. 

V  THE  NAVIER-STOKES  EQUATIONS  AND  CHARACTERISTIC  VARIABLES; 

We  now  begin  our  discussion  of  the  equations  of  gas  dynamics.  We 
will  neglect  viscosity  for  the  purposes  of  this  analysis.  We  will  assume 


that  the  flow  is  one-dimensional  and  subsonic  and  that  the  deviations 


from  free-stream  solutions  are  small.  This  will  allow  us  to  neglect 
second  order  terms. 

There  are  many  forms  of  this  equation,  but  the  one  most  suitable 
for  the  present  discussion  is 


It  ^ 


(6) 


where 


(y-  3)i 

wK  MU 


1 

(3-y)  u 

^  -  |(v-l)u^ 


and 


U  = 


or  in  terms  of  physical  variables 


3t  3x  ^ 


(7) 


where 


A  -  am 


and 


Here  we  make  the  key  assumption  that  deviations  from  the  free  stream 
are  going  to  be  sufficiently  small  that  we  can  treat  the  entries  in  the 
matrix  A  as  being  approximately  constant  (at  least  locally) .  Denote 
these  frozen  variables  by  a  0-subscript.  We  then  make  the  substitution 


and  when  this  is  substituted  into  (7)  we  obtain 


3Wi 

3t 

*^0  3x 

3W- 

8W 

4 

3t 

+ 

3W, 

3W 

3t 

+ 

<“o-=o>3r  -  “ 

Notice  now  how  this  breaks  down  into  two  separate  cases.  On  the  one  hand 
if  flow  is  supersonic  then  all  wave  motion  is  in  the  left  to  right  direc¬ 
tion.  In  this  case  all  analytic  boundary  conditions  ought  be  prescribed 
at  the  left  hand  side  and  only  numerical  boundary  conditions  prescribed 
at  the  right  hand  side. 

Since  the  substitution  (8)  is  equivalent  to 

P  •  Ki  +  (p^/2c^)(K2  +  K3) 

u  -  1/2  (K2  -  K3) 

P  •  (‘Si  +  ''3> 

it  follows  that  prescribing  all  physical  variables  at  the  inflow  and  pre¬ 
scribing  3p/3x  “  3u/3x  ■  3p/3x  ■  0  at  the  outflow  is  sound  in  terms  of 
analytical  and  numerical  requirements. 

However,  we  must  now  consider  the  case  of  subsonic  flow.  In  this 
case  the  situation  is  completely  different.  Here,  two  of  the  variables  W 
and  W2  go  left  to  right  with  velocities  u  and  u+c  respectively,  whereas 
one  of  the  variables  runs  right  to  left  with  velocity  c  -u  .  While 


the  variables  and  have  no  clear  physical  significance.  Yet  it 
is  only  by  considering  these  variables  that  the  full  wave  structure  of 
the  equations  (6)  or  (7)  can  be  understood.  Thus,  one  would  be  led  to 
predict,  for  small  deviations  from  free  stream  conditions,  that  the  best 
boundary  conditions  would  be  for  an  interval  (o,  L) 


(o.t)  - 


(o,t)  «  K2 


dW 


1  a,t) 


"^^2  (L,t)  -  0 
dx 


(9) 


^  (o,t)  -  0  W3  (L.t)  -  K3 

Note  the  curious  aspect  of  these  boundary  conditions.  In  order  to 
prescribe  the  numerical  values  and  K2,  we  need  to  know  accurately  all 
three  physical  variables  at  some  distance  to  the  left.  However,  only  the 
two  combinations  and  K2  are  prescribed.  This  can  be  summarized  by  saying 
that  while  we  have  used  all  three  pieces  of  information  upstream,  we  have 
done  so  in  such  a  way  that  one  degree  of  freedom  remains,  thus  allowing 
the  waves  in  to  exit  without  problems. 

On  the  basis  of  the  linearized  model,  various  other  combinations 
would  be  well-posed.  For  example,  it  is  possible  to  prescribe  in  terms 

of  either  or  K2  at  the  outflow  x  >  L.  Thus  at  the  outflow  one  may 
prescribe 


W2  (L.t)  -  F3  (t)  +  c^  Wj  (L.t)  +  C2  W2  (L.t) 

For  example  if  c^  ■  0.  C2  ■  1,  then  this  amounts  to  putting 

u  (L.t)  -  1/2  F3  (10) 

l.e.,  we  prescribe  velocity  at  the  outflow. 

Alternatively,  we  might  take  Cj^  -  0,  C2  -  -1  and  we  would  get 

p  (L,t)  -  ((p^c^)/2)F3  (11) 
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i.e.  we  prescribe  velocity  at  the  outflow.  Many  other  combinations  are 
possible,  but  as  remarked  in  section  IV,  all  these  will  cause  errors  in 
the  initial  data  to  be  reflected  back  into  the  medium  as  waves  running 


from  right  to  left.  For  example,  we  would  predict  that  an  error  in  W2 
would  be  reflected  back  as  an  error  in  if  we  use  boundary  condition 
(10).  As  we  shall  see,  this  is  exactly  what  happens. 

At  the  Inflow  end,  we  may  prescribe  and  in  terms  of  W^.  Thus 
the  following  boundary  conditions  are  well  posed; 

(o,t)  -  Fj  +  Cl  (o,t)  (12a) 

W2  (o,t)  -  F2 

For  example,  choosing  C2  *  +1  in  (Hh)  corresponds  to 
u  (o,t)  -  (1/2)  F2 


(i.e.  prescribing  u  at  the  Inflow)  and  C2  ■  -1  corresponds  to 
p  (o,t)  -  (p^Cj^/2)  F2 

(i.e.  prescribing  p) .  One  can  prescribe  the  combination  (u,p)  by  first 
choosing  C2  ■  1  (thereby  prescribing  u)  and  then  choosing  Cj^  «  p^/c^, 
thereby  prescribing  p  in  terms  of  a  given  F,  and  a  prescribed  u(o,t).  Sine 
(11a)  and  (11b)  reduce  to 


o  o 


p  -  vr  <Po/%  -  ‘^i^  ^^1“ '  ^1 

o  o 


about  the  only  condition  we  cannot  prescribe  is  u(o,t),  p(o,t),  since  there 


is  no  choice  of  c^,  C2  to  eliminate  p  from  these  equations. 

Again,  we  emphasize  that  each  of  these  boundary  conditions  is  reflect¬ 


ing,  i.e.  deviations  from  the  free  stream  in  the  initial  data  get  reflected 


back  as  waves  in  and  and  then  travel  back  downstream.  About  the 
worst  thing  that  can  be  done  is  to  prescribe  reflecting  boundary  condi¬ 
tions  at  the  inflow  x  «  0  and  the  outflow  x  «  L.  In  this  case  errors 
can  keep  being  reflected  up  and  down  the  region,  never  being  allowed  to 
exit.  This  prevents  convergence  to  a  steady  state  and  may  even  give  rise 
to  fictitious  periodic  oscillations. 

We  conclude  this  section  with  a  review  of  the  conclusions  on  bound¬ 
ary  conditions.  Based  on  the  linear  model,  boundary  conditions  (9)  seem 
optimal.  Any  other  prescription  of  the  physical  variables,  although  well 
posed,  causes  reflections  of  the  deviation  from  the  true  solution.  If  for 
example,  only  the  physical  variable  p^  is  known  at  the  outflow  then  it  is 
possible  to  prescribe  p  at  the  outflow,  in  such  a  manner  that  the  problem 
remains  well  posed.  We  emphasize  that  this  will  cause  errors  to  propagate 
upstream,  thereby  slowing  the  process  of  convergence  to  steady  state. 

This  may  not  be  too  bad,  so  long  as  the  upstream  boundary  conditions  are 
not  also  reflecting.  On  the  other  hand,  if  they  are,  then  convergence  to 
free  stream  may  never  occur. 


VI.  DISCUSSION  OF  NUMERICAL  RESULTS  -  NONLINEAR  COUPLING; 

The  one  dimensional  Navler-Stokes  equations  were  solved  with  an 
alternative  direction  explicit  MacCormack  scheme,  on  a  one-dimensional 
net  with  forty  grid  points.  The  code  was  an  exact  one-dimensional  version 
of  a  three-dimensional  code  which  had  proved  successful  in  many  supersonic 
studies  [6],  [ll].  There  were  two  questions  to  answer.  The  first  was, 
given  that  equation  (8)  is  correct  for  infinitesimally  small  deviations 
from  a  constant  free  steam,  how  correct  is  it  when  deviations  of  an 


Intermediate  (of  the  order  of  lOZ)  magnitude  are  there  instead?  It 
would  be  too  much  to  hope  that  the  variables  W^,  W^,  and  W^,  remain 
uncoupled,  but  we  should  be  able  to  get  an  idea  of  the  order  of  magni¬ 
tudes  of  the  coupling  involved.  The  second  problem  of  course,  is  to 
assess  the  Influence  of  the  various  types  of  boundary  conditions  commonly 
employed.  We  will  deal  with  the  latter  problem  in  section  VII. 

To  do  this,  we  considered  a  uniform  free  stream  situation,  with 

2 

pressure  equal  to  2000  Ibs/ft  ,  velocity  equal  to  548  ft/sec  and  density 

3 

equal  to  0.0023  slugs/ft  .  We  created  a  deviation  from  the  steady  state 
condition  in  a  variety  of  ways  as  In  [9]  ,  [lO],  and  then  watched  the 
progress  (or  lack  of  It)  to  a  steady  state.  A  variety  of  different  wrong 
initial  conditions  were  used.  One  type  was  to  impose  a  lOZ  deviation  in 
one  of  the  variables  K^,  K^,  at  the  points  {x^,  j-18,19k20,21,22}. 

We  would  then  watch  the  disturbance,  graphically  as  it  propagated  up  or 
down  stream.  Another  possibllty  was  to  put  in  uniformly  wrong  initial 
conditions  where  some  or  all  of  the  characteristic  variables  W^,  W^, 
are  perturbed  throughout  by  a  percentage  error  of  lOZ.  Each  plot  then 
showed  the  percentage  error,  with  the  different  curves  representing  the 
progress  of  time  as  one  ascends  the  plot.  The  curves  are  plotted  every 
twenty  five  time  steps  when  At  *  (.9) (Ax)/(u+c) .  We  also  point  out  the 
percentage  errors  in  W^,  W^,  and  at  the  end  of  the  run,  so  as  to  obtain 
information  on  relative  accuracy  and  speed  of  convergence  of  the  various 
methods. 

Figure  l*glves  the  results  of  an  experiment,  which  is  ideal  in  terms 
of  the  linear  theory.  An  initial  disturbance  in  the  left-running 

^Figures  are  located  at  end  of  report. 


cHafacteristlc  variable  is  given  (graphed  as  K^)  and  we  observe  deflec- 
tlons  In  the  variables  W^.  As  one  can  readily  see  from  the  pictures, 

the  disturbance  propagates  upstream  rapidly  until  convergence  is  reached 
(.12  agreement  with  physical  variables),  which  takes  place  within  130 
time  steps.  This  gives  us  six  curves.  Note  that  although  a  good  deal 
of  undershoot  and  overshoot  In  becomes  apparent,  there  Is  no  significant 
interaction  with  or  W^.  The  same  situation  appears  with  Initial  dis¬ 
turbances  In  W^,  the  slow-moving  right  running  wave.  From  this  experiment 
It  appears  that  deviations  In  either  or  will  not  effect  either  of 
the  other  two  variables.  However,  as  shown  In  Figures  2  and  3,  when  Initial 
disturbances  are  In  ^2’  ^  entirely  different  situation  exists.  Figure  2 
shows  a  uniform  Initial  disturbance  In  of  minus  ten  percent,  while 
and  are  left  undisturbed.  Initially  the  wave  In  W^,  propagates  rapidly 
out  of  the  medlvoa.  Indeed,  after  fifty  time  steps  It  Is  essentially  gone 
from  the  picture.  However,  this  does  not  happen  without  affecting  the  other 
two  variables.  Notice  how.  In  the  top  graph,  a  large  disturbance  Is  left 
in  after  fifty  time  steps  and  in  W^,  we  have  that  values  one  almost 
constant  at  minus  eight  percent.  However,  once  the  wave  has  made  Its 

exit,  the  other  two  variables  uncouple  and  resume  their  normal  wave  motion, 
and  the  error  can  be  seen  propagating  out  of  the  solution  In  the  usual  way 
convergence  Is  attained  within  three  hundred  Iterations.  These  particular 
results  illustrate  what  became  Increasingly  clear  throughout  the  series; 
perturbations  In  had  a  large  effect  on  and  whereas  If  was  not 
perturbed,  and  behaved  as  If  uncoupled.  Perturbations  in  and 
had  little  effect  on  H^.  No  qualitative  explanation  of  this  phenomenon 


Is  known  at  this  time. 


Figure  3  shows  the  same  phenomenon.  Here  an  error  of  -10%  is  made  in 

whereas  errors  of  +10%  are  made  in  and  Vy  Again,  we  see  that  unltl 
the  wave  exists,  there  are  massive  disturbances  in  the  wave  structure 
of  and  U^.  As  soon  as  exists,  (after  75  Interations)  the  regular 
wave  structlve  reasserts  Itself  and  errors  propagate  out  in  a  predictable 
wave-like  manner.  Again,  convergence  takes  approximately  three  hundred 
and  twenty  five  iterations.  It  seems  clear  that  this  is  optimal  given 
the  limited  wave  velocity,  so  we  can  ded\x:ethat  these  affects  are  due  to  the 
nonlinear  coupling.  Thus,  even  after  the  wave  exists  it  will  take  at 
least  a  mlximum  time  of  {L/(u-c),  L/u}  seconds  for  the  resulting  errors 
to  propagate  out  of  the  system. 


VII.  Discussion  OF  NUMERICAL  RESULTS  -  BOUNDARY  CONDITIONS; 

The  linear  "small  deflection"  theory  predicts  that  the  best  boundary 
conditions  would  be  the  prescription  of  the  charateristic  variables  at 
their  point  of  entry  with  some  form  of  (stable)  numerical  boundary  condi¬ 
tion  for  the  point  of  exit.  Here  are  two  such  schemes 

INFLOW  OUTFLOW 


W^  (o,t)  - 
(o,t)  -  K2 


(o,c) 


0 


Here  and  are  numbers  calculated  from  the  known  values  of 
u,  p,  p  at  the  inflow  and  is  calculated  from  the  known  values  of  u 
and  p  at  the  outflow.  Notice,  however,  the  one  "degree  of  freedom"  is 
left  at  the  Inflow  point.  This  allows  the  variables  to  adjust  but  in 
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compensating  ways.  The  boundary  conditions  In  the  code  are  usually 
In  terms  of  the  physical  variables  so  we  translate  (13)  to  physical 
variables. 


INFLOW 

Pi  -  [Kj  -  U2  -  (l/Po%)  P2] 

Uj^  -  1/2  [K2  +  ^2  -  (1/PpC^)  P2] 

^1  ‘  ^1  [K2  -  U2  +  a/p^c^)  P2] 


OUTFLOW 

“n  “  Cujj_^  +  a/Po^o^  Pn-1 

Pn  “  ^Po%/2)  [K3  +  Uj,_3  +  (1/p^c^) 

Pn  “  <Po^2c^)  [K3  +  Ujj.j]  -  [1/(2  C^2)]  +  Djj_i 

■^his  »et  of  boundary  conditions  is  predicted  to  work  well  in  the  Ifi.oar 
Siudies  of  one  equation,  occurring  in  [2]  and  [i].  We  sha]]  .  >!ll  .li.  se 
boundary  conditions  the  "no-change  characteristic  boundary  conditions". 
Another  possiblity,  suggested  by  one-D  analogues  in  [2],  is  the  following: 


INFLOW 

W^  (o.t)  -  K3 
W2  (o.t)  -  K2 


(u-c)^ 


0 


OUTFLOW 

5W, 
Uv  _ L 

-  0 

H 

(u+c)^ 

3t 

9x 

W3  (L,t)  -  K3 


(l^-<  ) 


where  derivatives  in  the  x-varlable  are  downwind  at  the  inflow  and  upwind 
at  the  outflow  and  forward  in  tine.  The  numbers  Kj^,  K2.  K3  are  prescribed 
as  before.  In  terms  of  the  physical  variables  these  translate  into 
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INFLOW 


-  1/2  [K^  +  u"  -d/p^c^)?; 

-  (p^c^/2)[K2+(1/p^c^)pJ-uJ(u^-c^) (At/Ax)  [uj-uj  +(l/PpC^) (pj+p")]] 
pj'*’^  -  Kj+(Pjj/2c^)Ck2+(1/p^c^)p"-uJ+(u^-c^)  (At/Ax)  [u2-u”  +(l/p^c^)  (pj-p”)  ]] 
OUTFLOW 

4^  -(l/2)[u"+^  pJ;-K3  +(Ax/Ax)(u^+c^)[u"_^-uJ+(1/p^c^)(p5_^+p")]] 
o  o 

Pn  “^Po‘=o/2)[K3+uJ  +(l/p^c^)p"+(At/Ax)(u^+c^)[uJ_3-u;;  ^PS-rPp^ 

pr^"^8^2c^)CK3+uJ  +(l/p^c^)pJ+(At/Ax)  (Vc^,)[u5_j-uJ  +(l/p^,c^) 

+  p"  -  (l/c^)p;  +  (At/Ax)u^Cp;_3-pJ+(l/c2)(p"  -  p“_^)] 

We  shall  call  these  the  "windward  difference  characterstic  boundary  condi¬ 
tions".  Note:  u  ,  c  can  be  different  values  at  inflow  and  outflow, 
o  o 

The  performance  of  the  code  with  either  of  these  was  analyzed  by  posing 
initial  conditions  in  which  there  was  a  disturbance  in  one  or  more  of  the 
characteristic  variables  either  locally  at  the  center  of  the  grid  or  uniform¬ 
ly  throughout  the  grid,  of  the  order  of  lOZ. 

Thus  figure  1  shows  what  happens  if  the  disturbance  in  only  in  the 
third  characteristics  variable  locally  using  the  windward  characteristic 
variables. 


Figure  3  shows  the  effect  of  a  plus  +10Z  error  in  the  initial  conditions 
W3  and  W3  and  a  -lOZ  error  in  W2.  (The  first  curve  from  the  bottom  is  the 
Initial  state  of  the  variable,  and  the  others  are  the  states  at  Intervals  of 
25  iterations).  In  figure  2,  we  have  an  initial  disturbance  in  W2  of  -10% 
with  no  initial  disturbance  in  W^^  or  W3.  The  pictures  look  essentially  the 
same  as  figure  1.  There  is  considerable  nonlinear  interaction  until  the 


W,  wave  exists,  and  then  uncoupled  wave  notion  to  the  right  In  the  first 
variable  (Wj^)  and  to  the  left  in  the  third  variable  (V^) .  There  are  no 
reflections  when  the  and  waves  exit  and  convergence  is  reached  in 
300  iterations. 

These  computations  were  made  using  the  windward  differencing  charac¬ 
teristic  boundary  conditions,  -Ithough  the  same  results  were  obtained  with 
the  no  change  characteristic  boundary  <  onditions. 

In  figure  4,  we  show  the  effect  of  a  local  disturbance  at  the  center 
of  the  grid  in  the  variable  with  the  second  set  of  B.C.'s  and  in  figure 
5  we  show  a  speeded  up  version  (every  50  Iterations)  of  the  same  distur¬ 
bance  with  the  first  set  of  B.C.'s.  The  third  and  fifth  graphs  on  figure  4 
are  almost  exactly  the  same  as  the  second  and  fourth  on  figure  5.  Figure  5 
shows  convergence  being  reached  in  300  iterations. 

We  conclude  that  either  of  the  first  two  sets  oi  boundary  conditions 
give  optimal  convergence  since  convergence  cannot  take  place  until  the  wave 
in  W2  exists  (very  quickly)  and  the  residual  (nonlinear)  effects  of  on 
can  exit  upstream.  If  they  can  do  this  without  any  reflections,  then 
the  convergence  is  essentially  optimal. 

Sometltres,  it  is  objected  that  in  a  wind  tunnel  experiment,  the  only 
variable  known  downstream  is  pressure  and  that  we  are  requiring  too  much 
information  in  prescribing  K^,  which  demands  a  knowledge  of  p  and  u  at 

the  outflow.  Suppose,  then  we  just  prescribe  at  the  outflow  using  the 

3W  dW 

otherwise  successful  conditions  of  _ 1  »  _ 2  ■  0  as  complementary  numerical 

9x  3x 

boundary  conditions.  Then  we  could  have,  for  example 
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ID/fLOh/ 

Wj^(o,t)  » 


OUTftOfc/ 


P  “  Pa 


W^Co.t)  =  K2 


sr  <’-■')  -  “ 


(o,t)  -  0 


=  (L,t)  -  0 


The  inflow  boundary  conditions  are  precisely  those  of  (12).  The  outflow 
boundary  conditions  (used  by  Steger  [l2])  are 


l>K  *  -  "s-l)  +  »I1-1 

“s  ■  -  fJ  *  "S-I 

The  predictions  of  section  VI  are  clear.  The  fact  that  p  is  pre¬ 
scribed  xeans  that  when  a  wave  in  comes  downstream,  it  exits  by  adjusting 
the  u  values  at  x  =  x„.  This  in  turn  causes  disturbances  in  W_*  -u+(l/p  c  )p 
which  cause  a  reflected  wave  upstream.  This  wave  can  be  seen  in  figure  6. 
Since  the  upstream  B.C.  is  non-reflecting,  this  means  that  the  left  running 
wave  will  exit  without  incidence.  When  compared  with  boundary  conditions 
(13)  or  (14)  it  is  obviously  less  desirable  because  of  the  magnitude  of  the 
reflection  in  W^.  However  it  does  converge  in  approximately  300  iterations, 
which  is  again  almost  optimal.  The  effect  of  these  large  oscillations  in 
more  complicated  geometries  may  prove  undesirable,  however. 

We  briefly  review  our  progress  so  far.  Two  sets  of  non-reflecting 
boundary  conditions  have  been  produced,  both  of  which  give  optimal  con¬ 
vergence  hut  which  rely  on  a  great  deal  of  information  at  both  ends.  The 
Infomation  however  is  used  in  such  a  way  as  to  allow  additional  degrees  of 
freedom  for  the  waves  to  exit  without  repeated  reflections.  One  reflecting 
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and  one  non -re flee ting  boundary  condition  can  be  combined  to  obtain 
almost  optimal  convergence  at  the  cost  of  some  large  left  running 
reflections,  whose  effect  in  more  complicated  geometries  remains  un¬ 
certain.  i-Thile  B.C.  13  was  expected  to  be  less  accurate  than  B.C.  14, 
little  evidence  for  this  has  been  uncovered,  except  at  the  boundaries. 

Time  dependent  periodic  flows  (e.g.  self  excited  oscillations)  may  prefer 
B.C.  14  however.  We  now  consider  some  of  the  other  boundary  condi¬ 
tions  which  have  been  tried  previously  in  the  literature. 

First  we  consider  the  case  of  reflecting  boundary  conditions.  Tliese 
occur  in  several  places  in  the  literature,  for  example  in  [lO]  and  [l2] 
and  have  been  discussed  in  section  VI,  As  we  have  seen,  these  arise  from 
prescribing  comblnatirtis  of  characteristic  variables  such  as  pressure  down- 
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St '•earn,  and  other  errf  ie.it ions  (perhaps 

to  density  and  ve1o''ity 

■J:  .jt  r.',>m)  , 

In  [iCj  Andy  .,nd  Strikwerda  considered 

(among  many  others)  the 

bcundary 

conditions 

INFLOW 

OUTFLOW 

u  »  u 

00 

dx  ° 

• 

T  »  T 

00 

0 

dx 

(16  ' 

dW 

=  0 
dx 

II 

8 

and  in  [l2]  Steger  uses 
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Figure  7  shows  the  effects  of  an  initial  local  disturbance  in  on 
both  of  these  sets  of  boundary  conditions,  (i6 )  on  the  left  and  (17)  on 
the  right.  Note  that  first  there  is  only  a  disturbance  in  the  bottom 
picture.  By  iteration  75  (fourth  curve  up  from  the  bottom)  we  can  see 
reflections  in  both  W2  and  although  the  deiration  is  a  little  harder 
to  see.  By  iteration  125  it  can  be  seen  that  the  wave  has  travelled 
downstream  and  is  reflected  back  upstream  in  W^.  These  reflections  con¬ 
tinued  (somewhat  smeared  out)  for  at  least  2000  iterations.  Perhaps  most 
startling  is  (at  least  at  the  beginning)  how  similar  they  are.  A  conclusion 
may  be  drawn  from  B.C.'s  (15;,  (16),  and  (17).  Prescribing  physical  values 
instead  of  characteristic  values  gives  rise  to  reflection.  If  the  reflec¬ 
tions  can  occur  only  at  one  end,  this  does  not  impede  convergence.  However, 
reflecting  conditions  at  both  ends  can  be  disastrous.  This  accounts  for 
the  "spurious  pressure  waves"  mentioned  by  Moretti  [?]. 

The  case  of  prescription  at  the  wrong  end  is  now  considered.  In  this 
case,  the  variables  u,  p,  p  are  prescribed  at  the  inflow  and  the  conditions 

dx  dx  dx 

are  prescribed  at  the  outflow.  In  this  case,  square  wave  disturbances 
which  affected  only  the  interiors  of  the  domain  exited  as  in  figure  1,  with 
some  minor  oscillations.  However,  when  a  uniformly  wrong  initial  condi¬ 
tion  was  imposed,  convergence  was  very  slow  with  large  oscillations,  and 
the  solution  converged  to  the  wrong  values,  with  errors  of  as  much  as  27%. 
The  converged  value  was  a  function  of  the  initial  condition  at  the  outflow 
end,  as  predicted  in  Gustafson  and  Krelss  [s].  The  resulting  graphs  are 


given  in  figure  8. 


boundaries.  This  r.ethod  is  •  entioned  in  [9]  as  giving  good  results  al¬ 
though  the  authors  caution  -sc  inst  it  on  the  grounds  of  siriall  "'S'  i  nations 
being  present.  In  fact  the  situation  is  much  more  serious.  If  initial 
waves  in  the  Interior  of  the  domain  are  used  with  the  initial  conditions 
correct  near  the  outflow,  then  the  solution  converges  rapidly  as  the 
travelling  vtaves  exit  without  reflections  and  with  minor  oscillations. 
However,  if  the  initial  data  is  uniformly  wrong  with  an  error  intially 
at  the  outflow  point,  then  and  converge  rapidly  but  accumulates 
huge  errors  of  the  order  of  70%.  Eventually  when  and  are  converged, 
the  correct  value  is  propagated  upwind  in  but  taking  large  amounts  of 
time  to  converge  because  of  the  large  errors  near  the  inflow’  point.  Indeed, 
as  bectiues  more  accurate  downstream  the  inaccuracies  become  much 
larger  (l'i0%)  upstream.  T-irtioularly  in  a  time  deper.duMU  problem  or  in 
a  problem  with  more  complicated  geometries,  this  could  he  truly  disastrous. 

It  points  to  another  fact:  if  a  new  boundary  condition  is  being  tested,  it 
is  not  sufficient  to  consider  initial  value  perturbations  from  free  stream 
which  are  non-zero  in  the  interior  only.  In  this  case,  we  might  have 
drawn  totally  wrong  conclusions  from  the  time  taken  for  convergence. 

In  [lO],  a  separate  non-reflecting  boundary  condition  is  proposed. 

This  boundary  condition  alters  the  value  of  K^,  increasing  it  if  the  computed 
value  of  p  is  less  than  p^,  and  decreasing  it  if  the  computed  value  of  p  is 
gi eater  than  p^.  Some  thought  shows  that  this  cannot  be  optimal.  Indeed, 
we  can  choose  a  variation  from  the  intial  conditions  in  which  p  is  less 
than  p  ,  but  because  us  is  smaller  than  u  ,  computed  V.  is  actually  larger 
than  In  this  case,  the  proposed  non-reflecting  boundary  condition 


of  [lO]  could  Introduce  more  errors  Into  the  system,  by  Increasing 
at  the  outflow  point  even  more. 

VIII .  RECOMMENDATIONS; 

We  have  completed  an  Initial  study  of  the  one-dimensional  Navler- 
Stokes  equations  and  their  far-fleld  boundary  conditions;  and  developed  . 
a  one-dimensional  code  to  provide  numerical  solutions.  We  have  used 
this  code  to  evaluate  the  Impact  of  a  variety  of  different  boundary  condi¬ 
tions  upon  the  successful  numerical  solution  of  these  problems.  Several 
outstanding  problems  remain  to  be  studied. 

First  one  should  examine  the  usefulness  of  the  two  recommended  boundary 
condition  sets  in  time  dependent  slttiatlons.  It  may  well  be  that  the 
difference  between  no-change  characteristic  boundary  conditions  and  windward- 
difference  characteristic  boundary  conditions  will  prove  to  be  highly  slgnlfl' 
cant  important  in  a  time  dependent  problem.  If  so,  this  would  be  signifi¬ 
cant  since  a  high  proportion  of  work  currently  done  at  the  Plight  Dynamics 
Laboratory  Computational  Aerodynamics  Group  is  of  this  type. 

Recommendation  for  future  efforts  is  to  Implement  these  characteristic 
boundary  conditions  in  two  and  three-dimensional  codes.  This  is  presently 
in  progress.  As  this  is  done,  there  is  no  doubt  that  further  refinements 
will  need  to  be  considered. 

In  the  discussion  of  our  results,  we  noted,  but  did  not  stress,  the 
following  problem.  Truly  accurate  time  dependent  solutions  only  exist  if 
the  C.F.L.  number  is  close  to  1.  The  time  step  is  limited  by  a  stability 
requirement  that  the  C.F.L.  number  be  less  than  1.  The  problem  with  this 
is  that  the  C.F.L.  number  is  calculated  on  the  basis  of  the  fastest  moving 


wave,  which  Is  the  fastest  to  leave  the  region.  After  a  relatively 
short  amount  of  time,  only  the  slow-moving  waves  are  left.  This  ought 
(at  least  for  some  cycle  of  Iterations)  allow  us  to  use  a  larger  time 
step.  This  has  not  yet  been  thoroughly  explored  and  we  recommend 
that  a  study  be  made  of  this  problem  This  could  result  In  substantial 
savings  of  expenditure  of  computer  time. 

More  work  needs  to  be  done  to  properly  understand  the  nonlinear 
actions  between  the  characteristic  waves.  This  problem  is  not  under¬ 
stood  at  this  time  and  a  proper  understanding  might  give  a  clue  to  correct 
posing  of  Intlal  conditions  In  such  a  way  as  to  minimize  start-up  shocks, 
as  observed  In  figures  4  and  5. 

Another  problem  to  be  Investigated  on  the  one-dimensional  code  Is 
whether  up-dating  the  coefficients  c^  u^  to  the  (n+l)  time  step 
values  significantly  accelerates  convergence  to  steady  state  or  improves 
accuracy.  Either  result  would  Improve  computational  efficiency  In  the 
future. 

Finally,  we  need  to  evaluate  how  appropriate  these  sets  of  boundary 
conditions  are  In  the  presence  of  stationary  shock  waves.  This  can  be 
done  first  on  a  one-dimensional  model  problem. 
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